close all
clear
%clc
addpath(genpath(pwd))
warning off
%% Section 1: User settings for the estimation
startYear     = 1961;       %Start year for estimation
endYear       = 2019;       %End year for estimation

lambdaSMM     = 1           %Weight in objective function to Shrinkage of based on Method of Moments (SMM)
SMMwithCSreg  = 1           %2 for incl. CS and risk-adj CS moments,
                            %1 for incl. CS moments,
                            %0 for not including CS moments

appOrder      = 3;          % Approximation order of DSGE model                            % 0 for using the mean in Perturbation to re-center the grid 

nx            = 6;          % #of states
ne            = 5;          % #of shocks

numCPUs       = feature('numcores'); % Number of CPUs

numOptimSteps = 0;          % Number of optimizations
optim         = 1;          % 1 Nelder-Mead.
                            % 2 Nelder-Mead with transformed parameter space
                            % 3 Own gradient based optimizer
                            % 4 CMAES
                            % 5 fmincon
cmaesPopSize  = 6*8;         % population size for CMAES
cmaesSigma    = 0.01;       % sigma in CMAES
cmaesSigmaMax = 0.05;        % sigmaMax in CMAES
MaxIter       = 50000;       % Max number of iterations for the optimizer
MaxEvals      = 50000;      % Max number of function evaluations for the optimizer
TolFun        = 1D-3;       % Function tolerance at the objective function
TolX          = 1D-3;       % Function tolerance for the coefficients
epsValue      = 1e-5;       % stepsize for computing gradient in own optimization routine
seOn          = 1;          % 1 for computing standard errors, else 0
ScalingMacro  = 0.2;        %Setting measurement errors for macro variables
ScalingYield  = 0.1;       %Setting measurement errors for yields

if numCPUs > 8
    % We use Grendel-S
    rmpath([pwd,'/Mex_myPC'])
    MEXon     = 1*0;          %1 For using MEX-files, else 0
else
    % We use My PC
    rmpath([pwd,'\Mex_GrendelS'])
    MEXon     = 1*0;
end

% Default setting
tau           = 400*0;      %Moments computed based on T*tau simulations, else closed-form solution used when tau = 0.
pruningOn     = 1;          %1 for using the pruned approximation, 0 for the unpruned approximation
nyMom         = 10;         %number of the first nyMom observables used for the SMM part
CSregress_m   = 4;          %2 for biannual implementation of CS regressions, 4 for annual etc.
inclSurveys   = 0;          %1 for including surveys, else 0

%% Section 2: Data and Model parameters
% Setting the calibrated parameters as in the sample ending in 2007Q4
structData = loadData_2020(startYear,endYear);
meanData   = nanmean(structData.data,1)';
stdData    = nanstd(structData.data,0,1)';
setCalibrateParams

CSdata     = CampbellSchillerRegression(structData.yieldsAll(1:end,:),CSregress_m);
%[CSdata.maturities;CSdata.CSbetta(2,:); CSdata.CSBetta_CI95Boot]
%[CSdata.maturities;CSdata.CSbetta(2,:);(CSdata.CSbetta(2,:)-1.96*CSdata.CSbetta_se(2,:));(CSdata.CSbetta(2,:)+1.96*CSdata.CSbetta_se(2,:))]
PredData    = predictiveRegression(structData.yieldsAll(1:end,:),CSregress_m);

% Compute empirical moments and weigthing matrix
dataSMM = momentsSMM(structData.data,structData.yieldsAll,nyMom,CSregress_m,SMMwithCSreg);
qLag = 3;
W    = getOptimalWeighting(qLag,dataSMM);

% The variables which we do not select for estimation are calibrated to the assigned values
calibrateParams.CRRA    = 10;

% Starting values
params.ALFA    = -20;
params.BETTA   = 0.995;
params.B       = 0.1;
params.PHI     = 0.075;
params.CHI     = 2.8;
params.XI      = 0.75;
params.KAPAw   = 0;
params.ETA     = 6;
params.PHIpai_1= 2;    %forward looking taylor - rule used
params.PHIc_1  = 0.0;
%params.PHIpai= 2;
%params.PHIc  = 0.0;
params.RHOz    = 0.98;
params.RHOd    = 0.99;
params.RHOn    = 0.95;
params.RHOp    = 0.99;
params.RHOa    = 0.99;
params.PAIss   = 1.0;
params.SIGMAmuz= 0.0004;
params.SIGMAd  = 0.03;
params.SIGMAn  = 0.03;
params.SIGMAp  = 0.001;
params.SIGMAa  = 0.01;


% With select = [2:10]*4/CSregress_m in SMM
calibrateParams        = rmfield(calibrateParams,'U0'); %use this to have U0 being active
if lambdaSMM  == 1 && SMMwithCSreg  == 1
    calibrateParams.ScalingYield = 0.083;   
    % Best specification with projection solution
    paramsValues = [-11.101942823607    0.988637503556    0.635426561457    0.089931047515    5.382449932453  ...
        0.380779218722    0.909884449017    8.801097086628    3.360386068027    1.034982531614  ...
        0.469664278984    0.984459501458    0.989512812916    0.995163270115    0.983275042959  ...
        1.023564607272    0.003720077294    0.128367989135    0.040220146959    0.000407800759  ...
        0.002139052407  ]';
    % Q =   -22.2404, done

    % Best specification with third order perturbation
    paramsValues = [-14.249698074180    0.948341430911    0.625449864341    0.168801100593   15.632860271997  ...
        0.395746291164    0.944856844456    9.773891262991    9.999728539786    2.939792646599  ...
        0.218131030872    0.998052769222    0.979902976539    0.987218069022    0.996179941277  ...
        1.023528786199    0.006088840104    0.052899786505    0.107965469851    0.000857510345  ...
        0.001483785409  ]';
    % Q=  -9.9486, done
    calibrateParams.PHIpai_1 = 9.999728539786;
    params = rmfield(params,'PHIpai_1');
    paramsValues = [-14.249698074180    0.948341430911    0.625449864341    0.168801100593   15.632860271997  ...
        0.395746291164    0.944856844456    9.773891262991                      2.939792646599  ...
        0.218131030872    0.998052769222    0.979902976539    0.987218069022    0.996179941277  ...
        1.023528786199    0.006088840104    0.052899786505    0.107965469851    0.000857510345  ...
        0.001483785409  ]';
end
    

selectParams = fieldnames(params);
params       = values2struct(paramsValues,selectParams);

%% Section 3: Setting setupFilter and setupProj
% The lower and upper bounds for the parameters
[lowerBounds,upperBounds,Insigma] = setLowerUpperBounds;
upperBounds.BETTA = 0.999;
% Constructing the struct setupFilter
constructSetupFilter
setupFilter.zlbON          = 1;

%% Section 3: Estimation
% tic
[Q0,errorMes,model0,~,outFilter0] = objectFuncQML(struc2values(params,setupFilter.selectParams),setupFilter);
Q0
% toc

if numOptimSteps > 0
    for i=1:numOptimSteps
        % We switch between the optimizer in optim and the Nelder Mead
        if mod(i,2) == 1
            setupFilter.optim = optim;
        else
            setupFilter.optim = 6;
        end
        outEstim = runOptimization(setupFilter,params);
        params = outEstim.paramsOpt;
        printToScreen(struct2array(params))
    end
else
    setupFilter.optim = 0;
    outEstim = runOptimization(setupFilter,params);
% Check closed form solution on long simulated sample 
%     setupFilter.numSim = 1D6;
%     shocks             = randn(5,setupFilter.numSim);
%     setupFilter.shocks = (shocks-mean(shocks,2))./repmat(std(shocks,0,2),1,setupFilter.numSim);
%     noise              = randn(4,setupFilter.numSim);
%     setupFilter.noise  = (noise-mean(noise,2))./repmat(std(noise,0,2),1,setupFilter.numSim);
%     noiseYieldsAll     = randn(40,setupFilter.numSim);
%     setupFilter.noiseYieldsAll  = (noiseYieldsAll-mean(noiseYieldsAll,2))./repmat(std(noiseYieldsAll,0,2),1,setupFilter.numSim);
%     outEstimLongSim = runOptimization(setupFilter,params);
%     [outEstimLongSim.outFilter.Q                        outEstim.outFilter.Q]
%     [outEstimLongSim.outFilter.momCondition'            outEstim.outFilter.momCondition']
%     [outEstimLongSim.outFilter.resMoments.meanYields    outEstim.outFilter.resMoments.meanYields]
%     [outEstimLongSim.outFilter.resMoments.stdYields     outEstim.outFilter.resMoments.stdYields]
%     [outEstimLongSim.outFilter.resMoments.meanForEstim  outEstim.outFilter.resMoments.meanForEstim]
%     [outEstimLongSim.outFilter.resMoments.stdForEstim   outEstim.outFilter.resMoments.stdForEstim]
end

%% Model diagnostics
outEstim.outFilter.Q
outDiagnostics = plotOutputFilter(outEstim);
plotTermPremia(outEstim)
plotCSloadings(CSdata,outEstim);
reportFit
CRRA = getCRRA(outEstim.model.params)
outEstim.newsDecomSim = inflVarianceRatio(outEstim,outEstim.outFilter.outSim.xAll_cu);
outEstim.newsDecom    = inflVarianceRatio(outEstim,outEstim.outFilter.xhat);
outEstim.newsDecomSim.summaryTable



%% Additional ex post estimation analysis
outEstim.UnconditionalMoments = unconMoments(outEstim.outFilter.ySimAsInData,outEstim.setupFilter);
outEstim.reducedForm          = doReducedFormAnalysis(outEstim);

pai0         = 0.05;
numSimTiming = 2000;
pathLength   = 5000;
outEstim.paiOpt = GetTimingPrem(outEstim,pai0,numSimTiming,pathLength);


% Standard naming of output file
filename = ['ModelPer',num2str(appOrder),'th_endYear',num2str(endYear),'_lambdaSMM',num2str(lambdaSMM),'_SMMwithCSreg',num2str(SMMwithCSreg),'_inclSurveys',num2str(inclSurveys)];
if seOn == 1
    filename = [filename,'_withSE'];
end
%save(filename,'outEstim','outDiagnostics','structData','CSdata')
